
% hidden causes and states
%==========================================================================
% x    - hidden states:
%   x.o(1) - oculomotor angle
%   x.o(2) - oculomotor velocity
%   x.x(1) - target angle - extrinsic coordinates
%
% v    - causal states: force on target
%
% g    - sensations:
%   g(1) - oculomotor angle (proprioception)
%   g(2) - oculomotor velocity
%   g(:) - visual input - intrinsic coordinates
%--------------------------------------------------------------------------
 
clear all
close all 

% Prepare Y data
%==========================================================================

fname = '../../../data/s3e2b4t4.asc';

edff = edffile(fname);
analysis = ET_DEM(edff);
[cells ea] = analysis.cellSignals();

Y = [ea cells];
g = analysis.g_function();

% Set-up
%==========================================================================

M(1).E.s = 1/2;                               % smoothness
M(1).E.n = 4;                                 % order of
M(1).E.d = 1;                                 % generalised motion
  
% angular frequency of target motion
%--------------------------------------------------------------------------
w  = 2*pi/(analysis.FpC/10); % Because of the downsampling
 
% sensory mappings with and without occlusion
%--------------------------------------------------------------------------
%g  = '[x.o; exp(-([-8:8]'' - x.x + x.o(1)).^2)*(x.x < 1/2)]';
%h  = '[x.o; exp(-([-8:8]'' - x.x + x.o(1)).^2)]';
 
 
% oculomotor latencies (sinusoidal movement)
%==========================================================================

% slow pursuit following with (second order) generative model
%--------------------------------------------------------------------------
x.o = [0.1;0.1];                                  % motor angle & velocity
x.x = 0.1;                                      % target location
%x = [0;0;0];
 
% level 1: Displacement dynamics and mapping to sensory/proprioception
%--------------------------------------------------------------------------
M(1).f = '[x.o(2); (v - x.o(1))/4 - x.o(2)/2; v - x.x]';

M(1).g = g;
M(1).x = x;                                   % hidden states
M(1).V = diag([exp(10),ones(1,17)*exp(10)]);                              % error precision
M(1).W = 1;exp(16);                              % error precision;
 
% level 2: With hidden (memory) states
%--------------------------------------------------------------------------
M(2).f  = '[x(2); -x(1)]*0.01';
M(2).g  = 'x(1)'; 
M(2).x  = [0.1;0.1];                             % hidden states
M(2).V  = exp(16);                             % error precision
M(2).W  = exp(1);                             % error precision
 
% level 3: Encoding frequency of memory states (U)
%--------------------------------------------------------------------------
% M(3).f = '[0]';
% M(3).g = 'x(1)';
% M(3).x = 1;
% %M(3).v = 8;
% M(3).W = exp(12);
% M(3).V = exp(16);
% 
% M(4).v = 0;
% M(4).V = exp(16);

DEM.M  = M;
DEM.Y  = Y';

%analysis.DEM = DEM;
%DEM.M(1).E.nD = 2;
%DEM.U = ones(1,numel(ea));
DEM.M(1).E.nE = 4;
DEM.M(1).E.nM = 4;
DEM = spm_DEM(DEM);
figure;
analysis.plotTrial();

%spm_figure('GetWin','Figure 2'); clf
%spm_dem_occlusion_movie(DEM)

%DEM = spm_DEM_generate(M,1*ones(1,100));

 